%%Solving for the transitional dynamics of a 3 country model
var rho_chi1   B1 gyy1 Y1 YT1 YNT1  C1 gt1 gtT1 gtNT1 Hr1 gz1  gp1 P1 PT1 PNT1 Pi1 PiT1 PiNT1 w1 Vinnov1  Z1 T1 TNT1 Y_L1 gy1 gc1 ;
var rho_chi2   B2 gyy2 Y2 YT2 YNT2  C2 gt2 gtT2 gtNT2 Hr2 gz2  gp2 P2 PT2 PNT2 Pi2 PiT2 PiNT2 w2 Vinnov2  Z2 T2 TNT2 Y_L2 gy2 gc2 ;
var rho_chi3   B3 gyy3 Y3 YT3 YNT3  C3 gt3 gtT3 gtNT3 Hr3 gz3  gp3 P3 PT3 PNT3 Pi3 PiT3 PiNT3 w3 Vinnov3  Z3 T3 TNT3 Y_L3 gy3 gc3 ;
var EX1 EX2 IM1 IM2 Yw r ;
var     chivar11 ga11  A11 pi_share11 pi_shareNT11   V11 Vinnovator11 rp11 ;
varexo distexo11 distexoNT11 chiexo11 ;
var     chivar12 ga12  A12 pi_share12 pi_shareNT12   V12 Vinnovator12 rp12 ;
varexo distexo12 distexoNT12 chiexo12 ;
var     chivar13 ga13  A13 pi_share13 pi_shareNT13   V13 Vinnovator13 rp13 ;
varexo distexo13 distexoNT13 chiexo13 ;
var     chivar21 ga21  A21 pi_share21 pi_shareNT21   V21 Vinnovator21 rp21 ;
varexo distexo21 distexoNT21 chiexo21 ;
var     chivar22 ga22  A22 pi_share22 pi_shareNT22   V22 Vinnovator22 rp22 ;
varexo distexo22 distexoNT22 chiexo22 ;
var     chivar23 ga23  A23 pi_share23 pi_shareNT23   V23 Vinnovator23 rp23 ;
varexo distexo23 distexoNT23 chiexo23 ;
var     chivar31 ga31  A31 pi_share31 pi_shareNT31   V31 Vinnovator31 rp31 ;
varexo distexo31 distexoNT31 chiexo31 ;
var     chivar32 ga32  A32 pi_share32 pi_shareNT32   V32 Vinnovator32 rp32 ;
varexo distexo32 distexoNT32 chiexo32 ;
var     chivar33 ga33  A33 pi_share33 pi_shareNT33   V33 Vinnovator33 rp33 ;
varexo distexo33 distexoNT33 chiexo33 ;
//@#if j != i
var   drp11 Ha11  J11 Jinnovator11  eps11;
//@#endif
//@#if j != i
var   drp12 Ha12  J12 Jinnovator12  eps12;
//@#endif
//@#if j != i
var   drp13 Ha13  J13 Jinnovator13  eps13;
//@#endif
//@#if j != i
var   drp21 Ha21  J21 Jinnovator21  eps21;
//@#endif
//@#if j != i
var   drp22 Ha22  J22 Jinnovator22  eps22;
//@#endif
//@#if j != i
var   drp23 Ha23  J23 Jinnovator23  eps23;
//@#endif
//@#if j != i
var   drp31 Ha31  J31 Jinnovator31  eps31;
//@#endif
//@#if j != i
var   drp32 Ha32  J32 Jinnovator32  eps32;
//@#endif
//@#if j != i
var   drp33 Ha33  J33 Jinnovator33  eps33;
//@#endif
// Model Parameters 
parameters  mbar betar sig betaa bet eta gam_barg;
parameters ipr1 ipr_rho1  Bbar1 lam1 L1  gam1 al1 LNT1 LT1 A1;
parameters ipr2 ipr_rho2  Bbar2 lam2 L2  gam2 al2 LNT2 LT2 A2;
parameters ipr3 ipr_rho3  Bbar3 lam3 L3  gam3 al3 LNT3 LT3 A3;
parameters   tau11 dist11 distNT11 epsbar11;
parameters   tau12 dist12 distNT12 epsbar12;
parameters   tau13 dist13 distNT13 epsbar13;
parameters   tau21 dist21 distNT21 epsbar21;
parameters   tau22 dist22 distNT22 epsbar22;
parameters   tau23 dist23 distNT23 epsbar23;
parameters   tau31 dist31 distNT31 epsbar31;
parameters   tau32 dist32 distNT32 epsbar32;
parameters   tau33 dist33 distNT33 epsbar33;
load Parameters.mat;
load Initial_values.mat;
//load Q.mat;
//load epsbar.mat;
mbar=1.25;
betar=0.5223;//0.4614;//0.3432;//0.2882;//0.4285;//0.4597;//0.3155;//692;//0.3931;//0.4622;//0.4450;//0.3883;//0.4661;//0.5904;//0.5918;//0.5497;//0.5301;//0.5424;//0.5562;//0.5519;//0.4825;//0.3;//0.36;//0.3316;//0.3673;//0.4999;
bet=0.985;
sig=5;
betaa=betar;//0.5;
eta=0.001;
gam_barg = 1;//0.3010;
//chi11 = chi(1,1);
///chi12 = chi(1,2);
//chi13 = chi(1,3);
//chi21 = chi(2,1);
//chi22 = chi(2,2);
//chi23 = chi(2,3);
//chi31 = chi(3,1);
//chi32 = chi(3,2);
//chi33 = chi(3,3);
ipr1 = ipr(1);
ipr2 = ipr(2);
ipr3 = ipr(3);
ipr_rho1 = ipr_rho(1);
ipr_rho2 = ipr_rho(2);
ipr_rho3 = ipr_rho(3);
lam1=lam(1);
lam2=lam(2);
lam3=lam(3);
Bbar1=0.;//Bbar(1);
Bbar2=0;//Bbar(2);
Bbar3=-0.;//Bbar(3);
A1=Q(1);
A2=Q(2);
A3=Q(3);
al1=alpha(1);
al2=alpha(2);
al3=alpha(3);
L1=L(1);
L2=L(2);
L3=L(3);
LT1=LT(1);
LT2=LT(2);
LT3=LT(3);
LNT1=LNT(1);
LNT2=LNT(2);
LNT3=LNT(3);
epsbar11=epsbar(1,1);
epsbar12=epsbar(1,2);
epsbar13=epsbar(1,3);
epsbar21=epsbar(2,1);
epsbar22=epsbar(2,2);
epsbar23=epsbar(2,3);
epsbar31=epsbar(3,1);
epsbar32=epsbar(3,2);
epsbar33=epsbar(3,3);
dist11=distT(1,1);
dist12=distT(1,2);
dist13=distT(1,3);
dist21=distT(2,1);
dist22=distT(2,2);
dist23=distT(2,3);
dist31=distT(3,1);
dist32=distT(3,2);
dist33=distT(3,3);
tau11=tau(1,1);
tau12=tau(1,2);
tau13=tau(1,3)/5;
tau21=tau(2,1);
tau22=tau(2,2);
tau23=tau(2,3);
tau31=tau(3,1);
tau32=tau(3,2);
tau33=tau(3,3);
distNT11=distNT(1,1);
distNT12=distNT(1,2);
distNT13=distNT(1,3);
distNT21=distNT(2,1);
distNT22=distNT(2,2);
distNT23=distNT(2,3);
distNT31=distNT(3,1);
distNT32=distNT(3,2);
distNT33=distNT(3,3);
gam1=lamNT(1);
gam2=lamNT(2);
gam3=lamNT(3);
model;
//Resource Constraint
exp(Y1)=exp(C1)+exp(Hr1)
+exp(Ha11)
+exp(Ha12)
+exp(Ha13)
;
exp(Y2)=exp(C2)+exp(Hr2)
+exp(Ha21)
+exp(Ha22)
+exp(Ha23)
;
exp(Y3)=exp(C3)+exp(Hr3)
+exp(Ha31)
+exp(Ha32)
+exp(Ha33)
;
//Euler Equation
rho_chi1 = 0.25*(ipr_rho1)^gam_barg;
1+eta*(B1*exp(Y1)-Bbar1)= bet*(1+gc1)^(-1)*exp(r);
rho_chi2 = 0.25*(ipr_rho2)^gam_barg;
1+eta*(B2*exp(Y2)-Bbar2)= bet*(1+gc2)^(-1)*exp(r);
rho_chi3 = 0.25*(ipr_rho3)^gam_barg;
1+eta*(B3*exp(Y3)-Bbar3)= bet*(1+gc3)^(-1)*exp(r);
B1*exp(Y1)+B2*exp(Y2)+B3*exp(Y3)=0;
//Prices
exp(PT1)^(1-sig)=
+((A1))^(sig-1)*(exp(T1(-1)))*(mbar*exp(w1)*dist11*(1+tau11*(1-distexo11)))^(1-sig)
+((A2))^(sig-1)*(exp(T2(-1)))*(mbar*exp(w2)*dist12*(1+tau12*(1-distexo12)))^(1-sig)
+((A3))^(sig-1)*(exp(T3(-1)))*(mbar*exp(w3)*dist13*(1+tau13*(1-distexo13)))^(1-sig)
;
exp(PT2)^(1-sig)=
+((A1))^(sig-1)*(exp(T1(-1)))*(mbar*exp(w1)*dist21*(1+tau21*(1-distexo21)))^(1-sig)
+((A2))^(sig-1)*(exp(T2(-1)))*(mbar*exp(w2)*dist22*(1+tau22*(1-distexo22)))^(1-sig)
+((A3))^(sig-1)*(exp(T3(-1)))*(mbar*exp(w3)*dist23*(1+tau23*(1-distexo23)))^(1-sig)
;
exp(PT3)^(1-sig)=
+((A1))^(sig-1)*(exp(T1(-1)))*(mbar*exp(w1)*dist31*(1+tau31*(1-distexo31)))^(1-sig)
+((A2))^(sig-1)*(exp(T2(-1)))*(mbar*exp(w2)*dist32*(1+tau32*(1-distexo32)))^(1-sig)
+((A3))^(sig-1)*(exp(T3(-1)))*(mbar*exp(w3)*dist33*(1+tau33*(1-distexo33)))^(1-sig)
;
exp(PNT1)^(1-sig)=
+((A1))^(sig-1)*exp(TNT1(-1))*(mbar*exp(w1)*distNT11*(1-distexoNT11))^(1-sig)
+((A2))^(sig-1)*exp(TNT2(-1))*(mbar*exp(w2)*distNT12*(1-distexoNT12))^(1-sig)
+((A3))^(sig-1)*exp(TNT3(-1))*(mbar*exp(w3)*distNT13*(1-distexoNT13))^(1-sig)
;
exp(PNT2)^(1-sig)=
+((A1))^(sig-1)*exp(TNT1(-1))*(mbar*exp(w1)*distNT21*(1-distexoNT21))^(1-sig)
+((A2))^(sig-1)*exp(TNT2(-1))*(mbar*exp(w2)*distNT22*(1-distexoNT22))^(1-sig)
+((A3))^(sig-1)*exp(TNT3(-1))*(mbar*exp(w3)*distNT23*(1-distexoNT23))^(1-sig)
;
exp(PNT3)^(1-sig)=
+((A1))^(sig-1)*exp(TNT1(-1))*(mbar*exp(w1)*distNT31*(1-distexoNT31))^(1-sig)
+((A2))^(sig-1)*exp(TNT2(-1))*(mbar*exp(w2)*distNT32*(1-distexoNT32))^(1-sig)
+((A3))^(sig-1)*exp(TNT3(-1))*(mbar*exp(w3)*distNT33*(1-distexoNT33))^(1-sig)
;
exp(P1)=(exp(PT1)/al1)^al1*(exp(PNT1)/(1-al1))^(1-al1);
exp(P2)=(exp(PT2)/al2)^al2*(exp(PNT2)/(1-al2))^(1-al2);
exp(P3)=(exp(PT3)/al3)^al3*(exp(PNT3)/(1-al3))^(1-al3);
//Trade shares
 exp(pi_share11)=((A1))^(sig-1)*exp(T1(-1))*(dist11*(1+tau11*(1-distexo11)))^(1-sig)*(exp(w1)*mbar)^(1-sig)/exp(PT1)^(1-sig);
 exp(pi_shareNT11)=((A1))^(sig-1)*exp(TNT1(-1))*(distNT11*(1-distexoNT11))^(1-sig)*(exp(w1)*mbar)^(1-sig)/exp(PNT1)^(1-sig);
 exp(pi_share12)=((A2))^(sig-1)*exp(T2(-1))*(dist12*(1+tau12*(1-distexo12)))^(1-sig)*(exp(w2)*mbar)^(1-sig)/exp(PT1)^(1-sig);
 exp(pi_shareNT12)=((A2))^(sig-1)*exp(TNT2(-1))*(distNT12*(1-distexoNT12))^(1-sig)*(exp(w2)*mbar)^(1-sig)/exp(PNT1)^(1-sig);
 exp(pi_share13)=((A3))^(sig-1)*exp(T3(-1))*(dist13*(1+tau13*(1-distexo13)))^(1-sig)*(exp(w3)*mbar)^(1-sig)/exp(PT1)^(1-sig);
 exp(pi_shareNT13)=((A3))^(sig-1)*exp(TNT3(-1))*(distNT13*(1-distexoNT13))^(1-sig)*(exp(w3)*mbar)^(1-sig)/exp(PNT1)^(1-sig);
 exp(pi_share21)=((A1))^(sig-1)*exp(T1(-1))*(dist21*(1+tau21*(1-distexo21)))^(1-sig)*(exp(w1)*mbar)^(1-sig)/exp(PT2)^(1-sig);
 exp(pi_shareNT21)=((A1))^(sig-1)*exp(TNT1(-1))*(distNT21*(1-distexoNT21))^(1-sig)*(exp(w1)*mbar)^(1-sig)/exp(PNT2)^(1-sig);
 exp(pi_share22)=((A2))^(sig-1)*exp(T2(-1))*(dist22*(1+tau22*(1-distexo22)))^(1-sig)*(exp(w2)*mbar)^(1-sig)/exp(PT2)^(1-sig);
 exp(pi_shareNT22)=((A2))^(sig-1)*exp(TNT2(-1))*(distNT22*(1-distexoNT22))^(1-sig)*(exp(w2)*mbar)^(1-sig)/exp(PNT2)^(1-sig);
 exp(pi_share23)=((A3))^(sig-1)*exp(T3(-1))*(dist23*(1+tau23*(1-distexo23)))^(1-sig)*(exp(w3)*mbar)^(1-sig)/exp(PT2)^(1-sig);
 exp(pi_shareNT23)=((A3))^(sig-1)*exp(TNT3(-1))*(distNT23*(1-distexoNT23))^(1-sig)*(exp(w3)*mbar)^(1-sig)/exp(PNT2)^(1-sig);
 exp(pi_share31)=((A1))^(sig-1)*exp(T1(-1))*(dist31*(1+tau31*(1-distexo31)))^(1-sig)*(exp(w1)*mbar)^(1-sig)/exp(PT3)^(1-sig);
 exp(pi_shareNT31)=((A1))^(sig-1)*exp(TNT1(-1))*(distNT31*(1-distexoNT31))^(1-sig)*(exp(w1)*mbar)^(1-sig)/exp(PNT3)^(1-sig);
 exp(pi_share32)=((A2))^(sig-1)*exp(T2(-1))*(dist32*(1+tau32*(1-distexo32)))^(1-sig)*(exp(w2)*mbar)^(1-sig)/exp(PT3)^(1-sig);
 exp(pi_shareNT32)=((A2))^(sig-1)*exp(TNT2(-1))*(distNT32*(1-distexoNT32))^(1-sig)*(exp(w2)*mbar)^(1-sig)/exp(PNT3)^(1-sig);
 exp(pi_share33)=((A3))^(sig-1)*exp(T3(-1))*(dist33*(1+tau33*(1-distexo33)))^(1-sig)*(exp(w3)*mbar)^(1-sig)/exp(PT3)^(1-sig);
 exp(pi_shareNT33)=((A3))^(sig-1)*exp(TNT3(-1))*(distNT33*(1-distexoNT33))^(1-sig)*(exp(w3)*mbar)^(1-sig)/exp(PNT3)^(1-sig);
//Labor market clearing condition
mbar*exp(w1)*LT1=
+exp(pi_share11)*exp(YT1)*exp(P1)/(1+tau11*(1-distexo11))
+exp(pi_share21)*exp(YT2)*exp(P2)/(1+tau21*(1-distexo21))
+exp(pi_share31)*exp(YT3)*exp(P3)/(1+tau31*(1-distexo31))
;
mbar*exp(w2)*LT2=
+exp(pi_share12)*exp(YT1)*exp(P1)/(1+tau12*(1-distexo12))
+exp(pi_share22)*exp(YT2)*exp(P2)/(1+tau22*(1-distexo22))
+exp(pi_share32)*exp(YT3)*exp(P3)/(1+tau32*(1-distexo32))
;
mbar*exp(w3)*LT3=
+exp(pi_share13)*exp(YT1)*exp(P1)/(1+tau13*(1-distexo13))
+exp(pi_share23)*exp(YT2)*exp(P2)/(1+tau23*(1-distexo23))
+exp(pi_share33)*exp(YT3)*exp(P3)/(1+tau33*(1-distexo33))
;
mbar*exp(w1)*LNT1=
+exp(pi_shareNT11)*exp(YNT1)*exp(P1)
+exp(pi_shareNT21)*exp(YNT2)*exp(P2)
+exp(pi_shareNT31)*exp(YNT3)*exp(P3)
;
mbar*exp(w2)*LNT2=
+exp(pi_shareNT12)*exp(YNT1)*exp(P1)
+exp(pi_shareNT22)*exp(YNT2)*exp(P2)
+exp(pi_shareNT32)*exp(YNT3)*exp(P3)
;
mbar*exp(w3)*LNT3=
+exp(pi_shareNT13)*exp(YNT1)*exp(P1)
+exp(pi_shareNT23)*exp(YNT2)*exp(P2)
+exp(pi_shareNT33)*exp(YNT3)*exp(P3)
;
// Output
exp(Y1)=exp(YT1);
exp(Y2)=exp(YT2);
exp(Y3)=exp(YT3);
exp(Yw)=exp(P1)*exp(Y1)+exp(P2)*exp(Y2)+exp(P3)*exp(Y3);
//Profits
exp(Pi1)= (1/(sig-1))*exp(w1)*L1; 
exp(PiT1)= (1/(sig-1))*exp(w1)*LT1; 
exp(PiNT1)= (1/(sig-1))*exp(w1)*LNT1; 
exp(Pi2)= (1/(sig-1))*exp(w2)*L2; 
exp(PiT2)= (1/(sig-1))*exp(w2)*LT2; 
exp(PiNT2)= (1/(sig-1))*exp(w2)*LNT2; 
exp(Pi3)= (1/(sig-1))*exp(w3)*L3; 
exp(PiT3)= (1/(sig-1))*exp(w3)*LT3; 
exp(PiNT3)= (1/(sig-1))*exp(w3)*LNT3; 
//Value of adopted technologies
exp(V11)=(1-exp(chivar11))*exp(PiT1)+(1/exp(r))*exp(V11(+1))*((1+gt3)^(1/(sig-1))/(1+(gt1)))*(exp(P1)/exp(P1(+1)));
exp(Vinnovator11)=exp(chivar11)*exp(PiT1)+(1/exp(r))*exp(Vinnovator11(+1))*((1+gt3)^(1/(sig-1))/(1+(gt1)))*(exp(P1)/exp(P1(+1)));
exp(chivar11) = (rho_chi1*exp(chiexo11));
exp(V12)=(1-exp(chivar12))*exp(PiT1)+(1/exp(r))*exp(V12(+1))*((1+gt3)^(1/(sig-1))/(1+(gt1)))*(exp(P1)/exp(P1(+1)));
exp(Vinnovator12)=exp(chivar12)*exp(PiT1)+(1/exp(r))*exp(Vinnovator12(+1))*((1+gt3)^(1/(sig-1))/(1+(gt1)))*(exp(P1)/exp(P1(+1)));
exp(chivar12) = (rho_chi1*exp(chiexo12));
exp(V13)=(1-exp(chivar13))*exp(PiT1)+(1/exp(r))*exp(V13(+1))*((1+gt3)^(1/(sig-1))/(1+(gt1)))*(exp(P1)/exp(P1(+1)));
exp(Vinnovator13)=exp(chivar13)*exp(PiT1)+(1/exp(r))*exp(Vinnovator13(+1))*((1+gt3)^(1/(sig-1))/(1+(gt1)))*(exp(P1)/exp(P1(+1)));
exp(chivar13) = (rho_chi1*exp(chiexo13));
exp(V21)=(1-exp(chivar21))*exp(PiT2)+(1/exp(r))*exp(V21(+1))*((1+gt3)^(1/(sig-1))/(1+(gt2)))*(exp(P2)/exp(P2(+1)));
exp(Vinnovator21)=exp(chivar21)*exp(PiT2)+(1/exp(r))*exp(Vinnovator21(+1))*((1+gt3)^(1/(sig-1))/(1+(gt2)))*(exp(P2)/exp(P2(+1)));
exp(chivar21) = (rho_chi2*exp(chiexo21));
exp(V22)=(1-exp(chivar22))*exp(PiT2)+(1/exp(r))*exp(V22(+1))*((1+gt3)^(1/(sig-1))/(1+(gt2)))*(exp(P2)/exp(P2(+1)));
exp(Vinnovator22)=exp(chivar22)*exp(PiT2)+(1/exp(r))*exp(Vinnovator22(+1))*((1+gt3)^(1/(sig-1))/(1+(gt2)))*(exp(P2)/exp(P2(+1)));
exp(chivar22) = (rho_chi2*exp(chiexo22));
exp(V23)=(1-exp(chivar23))*exp(PiT2)+(1/exp(r))*exp(V23(+1))*((1+gt3)^(1/(sig-1))/(1+(gt2)))*(exp(P2)/exp(P2(+1)));
exp(Vinnovator23)=exp(chivar23)*exp(PiT2)+(1/exp(r))*exp(Vinnovator23(+1))*((1+gt3)^(1/(sig-1))/(1+(gt2)))*(exp(P2)/exp(P2(+1)));
exp(chivar23) = (rho_chi2*exp(chiexo23));
exp(V31)=(1-exp(chivar31))*exp(PiT3)+(1/exp(r))*exp(V31(+1))*((1+gt3)^(1/(sig-1))/(1+(gt3)))*(exp(P3)/exp(P3(+1)));
exp(Vinnovator31)=exp(chivar31)*exp(PiT3)+(1/exp(r))*exp(Vinnovator31(+1))*((1+gt3)^(1/(sig-1))/(1+(gt3)))*(exp(P3)/exp(P3(+1)));
exp(chivar31) = (rho_chi3*exp(chiexo31));
exp(V32)=(1-exp(chivar32))*exp(PiT3)+(1/exp(r))*exp(V32(+1))*((1+gt3)^(1/(sig-1))/(1+(gt3)))*(exp(P3)/exp(P3(+1)));
exp(Vinnovator32)=exp(chivar32)*exp(PiT3)+(1/exp(r))*exp(Vinnovator32(+1))*((1+gt3)^(1/(sig-1))/(1+(gt3)))*(exp(P3)/exp(P3(+1)));
exp(chivar32) = (rho_chi3*exp(chiexo32));
exp(V33)=(1-exp(chivar33))*exp(PiT3)+(1/exp(r))*exp(V33(+1))*((1+gt3)^(1/(sig-1))/(1+(gt3)))*(exp(P3)/exp(P3(+1)));
exp(Vinnovator33)=exp(chivar33)*exp(PiT3)+(1/exp(r))*exp(Vinnovator33(+1))*((1+gt3)^(1/(sig-1))/(1+(gt3)))*(exp(P3)/exp(P3(+1)));
exp(chivar33) = (rho_chi3*exp(chiexo33));
//Probability of adoption
exp(eps11)=epsbar11*(exp(P1)*exp(Ha11)/exp(Yw))^betaa;
exp(eps12)=epsbar12*(exp(P1)*exp(Ha12)/exp(Yw))^betaa;
exp(eps13)=epsbar13*(exp(P1)*exp(Ha13)/exp(Yw))^betaa;
exp(eps21)=epsbar21*(exp(P2)*exp(Ha21)/exp(Yw))^betaa;
exp(eps22)=epsbar22*(exp(P2)*exp(Ha22)/exp(Yw))^betaa;
exp(eps23)=epsbar23*(exp(P2)*exp(Ha23)/exp(Yw))^betaa;
exp(eps31)=epsbar31*(exp(P3)*exp(Ha31)/exp(Yw))^betaa;
exp(eps32)=epsbar32*(exp(P3)*exp(Ha32)/exp(Yw))^betaa;
exp(eps33)=epsbar33*(exp(P3)*exp(Ha33)/exp(Yw))^betaa;
//Value of an innovation
exp(Vinnov1)=
   +exp(Jinnovator11)*exp(T1(-1))/exp(T1(-1))
   +exp(Jinnovator21)*exp(T1(-1))/exp(T2(-1))
   +exp(Jinnovator31)*exp(T1(-1))/exp(T3(-1))
;
exp(Vinnov2)=
   +exp(Jinnovator12)*exp(T2(-1))/exp(T1(-1))
   +exp(Jinnovator22)*exp(T2(-1))/exp(T2(-1))
   +exp(Jinnovator32)*exp(T2(-1))/exp(T3(-1))
;
exp(Vinnov3)=
   +exp(Jinnovator13)*exp(T3(-1))/exp(T1(-1))
   +exp(Jinnovator23)*exp(T3(-1))/exp(T2(-1))
   +exp(Jinnovator33)*exp(T3(-1))/exp(T3(-1))
;
//FOC innovation    
exp(Hr1)=exp(P1)^(-1)*(exp(Vinnov1)*betar*lam1*(exp(Yw))^(-betar))^(1/(1-betar));
exp(Hr2)=exp(P2)^(-1)*(exp(Vinnov2)*betar*lam2*(exp(Yw))^(-betar))^(1/(1-betar));
exp(Hr3)=exp(P3)^(-1)*(exp(Vinnov3)*betar*lam3*(exp(Yw))^(-betar))^(1/(1-betar));
//Law of motion of innovation
gz1=lam1*(exp(P1)*exp(Hr1)/exp(Yw))^betar*(exp(T1(-1))/exp(Z1(-1)));
gz2=lam2*(exp(P2)*exp(Hr2)/exp(Yw))^betar*(exp(T2(-1))/exp(Z2(-1)));
gz3=lam3*(exp(P3)*exp(Hr3)/exp(Yw))^betar*(exp(T3(-1))/exp(Z3(-1)));
//Law of motion of adoption
ga11=exp(eps11)*(exp(Z1(-1))*(1+gz1)/exp(A11(-1))-1);
ga12=exp(eps12)*(exp(Z2(-1))*(1+gz2)/exp(A12(-1))-1);
ga13=exp(eps13)*(exp(Z3(-1))*(1+gz3)/exp(A13(-1))-1);
ga21=exp(eps21)*(exp(Z1(-1))*(1+gz1)/exp(A21(-1))-1);
ga22=exp(eps22)*(exp(Z2(-1))*(1+gz2)/exp(A22(-1))-1);
ga23=exp(eps23)*(exp(Z3(-1))*(1+gz3)/exp(A23(-1))-1);
ga31=exp(eps31)*(exp(Z1(-1))*(1+gz1)/exp(A31(-1))-1);
ga32=exp(eps32)*(exp(Z2(-1))*(1+gz2)/exp(A32(-1))-1);
ga33=exp(eps33)*(exp(Z3(-1))*(1+gz3)/exp(A33(-1))-1);
//Value of an unadopted technology
exp(J11)=-exp(P1)*exp(Ha11)*(exp(T1)/exp(A11(-1)))*(exp(eps11)/(ga11))+(1/exp(r))*(exp(eps11)*exp(V11(+1))+(1-exp(eps11))*exp(J11(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt1)))*(exp(P1)/exp(P1(+1)));
exp(Jinnovator11)=(1/exp(r))*(exp(eps11)*exp(Vinnovator11(+1))+(1-exp(eps11))*exp(Jinnovator11(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt1)))*(exp(P1)/exp(P1(+1)));
exp(J12)=-exp(P1)*exp(Ha12)*(exp(T1)/exp(A12(-1)))*(exp(eps12)/(ga12))+(1/exp(r))*(exp(eps12)*exp(V12(+1))+(1-exp(eps12))*exp(J12(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt1)))*(exp(P1)/exp(P1(+1)));
exp(Jinnovator12)=(1/exp(r))*(exp(eps12)*exp(Vinnovator12(+1))+(1-exp(eps12))*exp(Jinnovator12(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt1)))*(exp(P1)/exp(P1(+1)));
exp(J13)=-exp(P1)*exp(Ha13)*(exp(T1)/exp(A13(-1)))*(exp(eps13)/(ga13))+(1/exp(r))*(exp(eps13)*exp(V13(+1))+(1-exp(eps13))*exp(J13(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt1)))*(exp(P1)/exp(P1(+1)));
exp(Jinnovator13)=(1/exp(r))*(exp(eps13)*exp(Vinnovator13(+1))+(1-exp(eps13))*exp(Jinnovator13(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt1)))*(exp(P1)/exp(P1(+1)));
exp(J21)=-exp(P2)*exp(Ha21)*(exp(T2)/exp(A21(-1)))*(exp(eps21)/(ga21))+(1/exp(r))*(exp(eps21)*exp(V21(+1))+(1-exp(eps21))*exp(J21(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt2)))*(exp(P2)/exp(P2(+1)));
exp(Jinnovator21)=(1/exp(r))*(exp(eps21)*exp(Vinnovator21(+1))+(1-exp(eps21))*exp(Jinnovator21(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt2)))*(exp(P2)/exp(P2(+1)));
exp(J22)=-exp(P2)*exp(Ha22)*(exp(T2)/exp(A22(-1)))*(exp(eps22)/(ga22))+(1/exp(r))*(exp(eps22)*exp(V22(+1))+(1-exp(eps22))*exp(J22(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt2)))*(exp(P2)/exp(P2(+1)));
exp(Jinnovator22)=(1/exp(r))*(exp(eps22)*exp(Vinnovator22(+1))+(1-exp(eps22))*exp(Jinnovator22(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt2)))*(exp(P2)/exp(P2(+1)));
exp(J23)=-exp(P2)*exp(Ha23)*(exp(T2)/exp(A23(-1)))*(exp(eps23)/(ga23))+(1/exp(r))*(exp(eps23)*exp(V23(+1))+(1-exp(eps23))*exp(J23(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt2)))*(exp(P2)/exp(P2(+1)));
exp(Jinnovator23)=(1/exp(r))*(exp(eps23)*exp(Vinnovator23(+1))+(1-exp(eps23))*exp(Jinnovator23(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt2)))*(exp(P2)/exp(P2(+1)));
exp(J31)=-exp(P3)*exp(Ha31)*(exp(T3)/exp(A31(-1)))*(exp(eps31)/(ga31))+(1/exp(r))*(exp(eps31)*exp(V31(+1))+(1-exp(eps31))*exp(J31(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt3)))*(exp(P3)/exp(P3(+1)));
exp(Jinnovator31)=(1/exp(r))*(exp(eps31)*exp(Vinnovator31(+1))+(1-exp(eps31))*exp(Jinnovator31(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt3)))*(exp(P3)/exp(P3(+1)));
exp(J32)=-exp(P3)*exp(Ha32)*(exp(T3)/exp(A32(-1)))*(exp(eps32)/(ga32))+(1/exp(r))*(exp(eps32)*exp(V32(+1))+(1-exp(eps32))*exp(J32(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt3)))*(exp(P3)/exp(P3(+1)));
exp(Jinnovator32)=(1/exp(r))*(exp(eps32)*exp(Vinnovator32(+1))+(1-exp(eps32))*exp(Jinnovator32(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt3)))*(exp(P3)/exp(P3(+1)));
exp(J33)=-exp(P3)*exp(Ha33)*(exp(T3)/exp(A33(-1)))*(exp(eps33)/(ga33))+(1/exp(r))*(exp(eps33)*exp(V33(+1))+(1-exp(eps33))*exp(J33(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt3)))*(exp(P3)/exp(P3(+1)));
exp(Jinnovator33)=(1/exp(r))*(exp(eps33)*exp(Vinnovator33(+1))+(1-exp(eps33))*exp(Jinnovator33(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt3)))*(exp(P3)/exp(P3(+1)));
//FOC adoption
exp(P1)*exp(Ha11)*(exp(T1(-1))/exp(A11(-1)))*(exp(eps11)/(ga11))=betaa*(1/exp(r))*exp(eps11)*(exp(V11(+1))-exp(J11(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt1)))*(exp(P1)/exp(P1(+1)));
exp(P1)*exp(Ha12)*(exp(T1(-1))/exp(A12(-1)))*(exp(eps12)/(ga12))=betaa*(1/exp(r))*exp(eps12)*(exp(V12(+1))-exp(J12(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt1)))*(exp(P1)/exp(P1(+1)));
exp(P1)*exp(Ha13)*(exp(T1(-1))/exp(A13(-1)))*(exp(eps13)/(ga13))=betaa*(1/exp(r))*exp(eps13)*(exp(V13(+1))-exp(J13(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt1)))*(exp(P1)/exp(P1(+1)));
exp(P2)*exp(Ha21)*(exp(T2(-1))/exp(A21(-1)))*(exp(eps21)/(ga21))=betaa*(1/exp(r))*exp(eps21)*(exp(V21(+1))-exp(J21(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt2)))*(exp(P2)/exp(P2(+1)));
exp(P2)*exp(Ha22)*(exp(T2(-1))/exp(A22(-1)))*(exp(eps22)/(ga22))=betaa*(1/exp(r))*exp(eps22)*(exp(V22(+1))-exp(J22(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt2)))*(exp(P2)/exp(P2(+1)));
exp(P2)*exp(Ha23)*(exp(T2(-1))/exp(A23(-1)))*(exp(eps23)/(ga23))=betaa*(1/exp(r))*exp(eps23)*(exp(V23(+1))-exp(J23(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt2)))*(exp(P2)/exp(P2(+1)));
exp(P3)*exp(Ha31)*(exp(T3(-1))/exp(A31(-1)))*(exp(eps31)/(ga31))=betaa*(1/exp(r))*exp(eps31)*(exp(V31(+1))-exp(J31(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt3)))*(exp(P3)/exp(P3(+1)));
exp(P3)*exp(Ha32)*(exp(T3(-1))/exp(A32(-1)))*(exp(eps32)/(ga32))=betaa*(1/exp(r))*exp(eps32)*(exp(V32(+1))-exp(J32(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt3)))*(exp(P3)/exp(P3(+1)));
exp(P3)*exp(Ha33)*(exp(T3(-1))/exp(A33(-1)))*(exp(eps33)/(ga33))=betaa*(1/exp(r))*exp(eps33)*(exp(V33(+1))-exp(J33(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt3)))*(exp(P3)/exp(P3(+1)));
//Total stock of knowledge
exp(T1)=
   +exp(A11)
   +exp(A12)
   +exp(A13)
;
exp(T2)=
   +exp(A21)
   +exp(A22)
   +exp(A23)
;
exp(T3)=
   +exp(A31)
   +exp(A32)
   +exp(A33)
;
//Growth rates
(gp1)=(exp(P1(+1))/exp(P1)-1)+((1/(1-sig))*gt3);
exp(TNT1)=gam1*exp(T1);
(gp2)=(exp(P2(+1))/exp(P2)-1)+((1/(1-sig))*gt3);
exp(TNT2)=gam2*exp(T2);
(gp3)=(exp(P3(+1))/exp(P3)-1)+((1/(1-sig))*gt3);
exp(TNT3)=gam3*exp(T3);
(gz1)=(exp(Z1)/exp(Z1(-1))-1)+gt3;
(gz2)=(exp(Z2)/exp(Z2(-1))-1)+gt3;
(gz3)=(exp(Z3)/exp(Z3(-1))-1)+gt3;
(gt1)=(exp(T1)/exp(T1(-1))-1)+gt3;
(gt2)=(exp(T2)/exp(T2(-1))-1)+gt3;
(gt3)=exp(A31(-1))*(ga31)+exp(A32(-1))*(ga32)+exp(A33(-1))*(ga33);
(gtT1)=gt1;
(gtNT1)=gt1;
(gtT2)=gt2;
(gtNT2)=gt2;
(gtT3)=gt3;
(gtNT3)=gt3;
(gy1)=((exp(Y_L1(+1)))/(exp(Y_L1))-1)+(1+gt3)^(1/(sig-1))-1;
(gyy1)=((exp(Y1(+1)))/(exp(Y1))-1)+(1+gt3)^(1/(sig-1))-1;
(gy2)=((exp(Y_L2(+1)))/(exp(Y_L2))-1)+(1+gt3)^(1/(sig-1))-1;
(gyy2)=((exp(Y2(+1)))/(exp(Y2))-1)+(1+gt3)^(1/(sig-1))-1;
(gy3)=((exp(Y_L3(+1)))/(exp(Y_L3))-1)+(1+gt3)^(1/(sig-1))-1;
(gyy3)=((exp(Y3(+1)))/(exp(Y3))-1)+(1+gt3)^(1/(sig-1))-1;
(drp11)=(exp(rp11(+1))/exp(rp11)-1)+(1+gt3)^(1/(sig-1))-1;
(drp12)=(exp(rp12(+1))/exp(rp12)-1)+(1+gt3)^(1/(sig-1))-1;
(drp13)=(exp(rp13(+1))/exp(rp13)-1)+(1+gt3)^(1/(sig-1))-1;
(drp21)=(exp(rp21(+1))/exp(rp21)-1)+(1+gt3)^(1/(sig-1))-1;
(drp22)=(exp(rp22(+1))/exp(rp22)-1)+(1+gt3)^(1/(sig-1))-1;
(drp23)=(exp(rp23(+1))/exp(rp23)-1)+(1+gt3)^(1/(sig-1))-1;
(drp31)=(exp(rp31(+1))/exp(rp31)-1)+(1+gt3)^(1/(sig-1))-1;
(drp32)=(exp(rp32(+1))/exp(rp32)-1)+(1+gt3)^(1/(sig-1))-1;
(drp33)=(exp(rp33(+1))/exp(rp33)-1)+(1+gt3)^(1/(sig-1))-1;
(gc1)=(exp(C1(+1))/exp(C1)-1)+(1+gt3)^(1/(sig-1))-1;
(gc2)=(exp(C2(+1))/exp(C2)-1)+(1+gt3)^(1/(sig-1))-1;
(gc3)=(exp(C3(+1))/exp(C3)-1)+(1+gt3)^(1/(sig-1))-1;
(ga11)=(exp(A11)/exp(A11(-1))-1)+gt3;
(ga12)=(exp(A12)/exp(A12(-1))-1)+gt3;
(ga13)=(exp(A13)/exp(A13(-1))-1)+gt3;
(ga21)=(exp(A21)/exp(A21(-1))-1)+gt3;
(ga22)=(exp(A22)/exp(A22(-1))-1)+gt3;
(ga23)=(exp(A23)/exp(A23(-1))-1)+gt3;
(ga31)=(exp(A31)/exp(A31(-1))-1)+gt3;
(ga32)=(exp(A32)/exp(A32(-1))-1)+gt3;
(ga33)=(exp(A33)/exp(A33(-1))-1)+gt3;
//Royalties
exp(rp11)=(exp(chivar11))*exp(PiT1)*exp(A11(-1))/exp(T1(-1));
exp(rp12)=(exp(chivar12))*exp(PiT1)*exp(A12(-1))/exp(T1(-1));
exp(rp13)=(exp(chivar13))*exp(PiT1)*exp(A13(-1))/exp(T1(-1));
exp(rp21)=(exp(chivar21))*exp(PiT2)*exp(A21(-1))/exp(T2(-1));
exp(rp22)=(exp(chivar22))*exp(PiT2)*exp(A22(-1))/exp(T2(-1));
exp(rp23)=(exp(chivar23))*exp(PiT2)*exp(A23(-1))/exp(T2(-1));
exp(rp31)=(exp(chivar31))*exp(PiT3)*exp(A31(-1))/exp(T3(-1));
exp(rp32)=(exp(chivar32))*exp(PiT3)*exp(A32(-1))/exp(T3(-1));
exp(rp33)=(exp(chivar33))*exp(PiT3)*exp(A33(-1))/exp(T3(-1));
//Trade balance equation
exp(EX1)=
+exp(pi_share21)*exp(P2)*exp(YT2)/(1+tau21*(1-distexo21))+exp(rp21)
+exp(pi_share31)*exp(P3)*exp(YT3)/(1+tau31*(1-distexo31))+exp(rp31)
;
exp(EX2)=
+exp(pi_share12)*exp(P1)*exp(YT1)/(1+tau12*(1-distexo12))+exp(rp12)
+exp(pi_share32)*exp(P3)*exp(YT3)/(1+tau32*(1-distexo32))+exp(rp32)
;
exp(IM1)=
+exp(pi_share12)*exp(P1)*exp(YT1)/(1+tau12*(1-distexo12))+exp(rp12)
+exp(pi_share13)*exp(P1)*exp(YT1)/(1+tau13*(1-distexo13))+exp(rp13)
;
exp(IM2)=
+exp(pi_share21)*exp(P2)*exp(YT2)/(1+tau21*(1-distexo21))+exp(rp21)
+exp(pi_share23)*exp(P2)*exp(YT2)/(1+tau23*(1-distexo23))+exp(rp23)
;
exp(IM1)=exp(EX1)+B1(-1)*exp(r(-1))*exp(Y1(-1))*(1/exp(gy1))-B1*exp(Y1);
exp(IM2)=exp(EX2)+B2(-1)*exp(r(-1))*exp(Y2(-1))*(1/exp(gy2))-B2*exp(Y2);
(w3)=(0);
//Productivity
  exp(Y_L1)=exp(Y1);
  exp(Y_L2)=exp(Y2);
  exp(Y_L3)=exp(Y3);
 
end;
initval;
B1=0;
B2=0;
B3=0;
distexo11=0;
distexo12=0;
distexo13=0;
distexo21=0;
distexo22=0;
distexo23=0;
distexo31=0;
distexo32=0;
distexo33=0;
distexoNT11=0;
distexoNT12=0;
distexoNT13=0;
distexoNT21=0;
distexoNT22=0;
distexoNT23=0;
distexoNT31=0;
distexoNT32=0;
distexoNT33=0;
chiexo11=0.;
chiexo12=0.0;
chiexo13=0.0;
chiexo21=0.;
chiexo22=0.0;
chiexo23=0.0;
chiexo31=0.;
chiexo32=0.0;
chiexo33=0.0;
r=log((1+ss_g/(sig-1))/(1+ss_g/(1-sig))/beta);
C1=log(ss_c(1));
C2=log(ss_c(2));
C3=log(ss_c(3));
Y1=log(ss_y(1));
Y2=log(ss_y(2));
Y3=log(ss_y(3));
YT1=log(ss_yT(1));
YT2=log(ss_yT(2));
YT3=log(ss_yT(3));
YNT1=log(ss_yNT(1));
YNT2=log(ss_yNT(2));
YNT3=log(ss_yNT(3));
T1=log(ss_TT(1)/ss_TT(3));
T2=log(ss_TT(2)/ss_TT(3));
T3=log(ss_TT(3)/ss_TT(3));
TNT1=log(ss_TNT(1)/ss_TNT(3));
TNT2=log(ss_TNT(2)/ss_TNT(3));
TNT3=log(ss_TNT(3)/ss_TNT(3));
P1=log(ss_p(1));
P2=log(ss_p(2));
P3=log(ss_p(3));
PT1=log(ss_pT(1));
PT2=log(ss_pT(2));
PT3=log(ss_pT(3));
PNT1=log(ss_pNT(1));
PNT2=log(ss_pNT(2));
PNT3=log(ss_pNT(3));
EX1=log(ss_EX(1));
EX2=log(ss_EX(2));
IM1=log(ss_IM(1));
IM2=log(ss_IM(2));
Pi1=log(ss_pi(1));
Pi2=log(ss_pi(2));
Pi3=log(ss_pi(3));
PiT1=log(ss_piT(1));
PiT2=log(ss_piT(2));
PiT3=log(ss_piT(3));
PiNT1=log(ss_piNT(1));
PiNT2=log(ss_piNT(2));
PiNT3=log(ss_piNT(3));
w1=log(ss_w(1));
w2=log(ss_w(2));
w3=log(ss_w(3));
Vinnov1=log(ss_vinnov(1));
Vinnov2=log(ss_vinnov(2));
Vinnov3=log(ss_vinnov(3));
chivar11 = log(ss_chi(1,1));
chivar12 = log(ss_chi(1,2));
chivar13 = log(ss_chi(1,3));
chivar21 = log(ss_chi(2,1));
chivar22 = log(ss_chi(2,2));
chivar23 = log(ss_chi(2,3));
chivar31 = log(ss_chi(3,1));
chivar32 = log(ss_chi(3,2));
chivar33 = log(ss_chi(3,3));
Z1=log(ss_Z(1));
Z2=log(ss_Z(2));
Z3=log(ss_Z(3));
Y_L1=log(ss_YL(1));
Y_L2=log(ss_YL(2));
Y_L3=log(ss_YL(3));
Hr1=log(ss_hr(1));
Hr2=log(ss_hr(2));
Hr3=log(ss_hr(3));
gz1=ss_g;
gz2=ss_g;
gz3=ss_g;
gt1=ss_g;
gt2=ss_g;
gt3=ss_g;
gtT1=ss_g;
gtT2=ss_g;
gtT3=ss_g;
gtNT1=ss_g;
gtNT2=ss_g;
gtNT3=ss_g;
gy1=ss_g/(sig-1);
gy2=ss_g/(sig-1);
gy3=ss_g/(sig-1);
gc1=ss_g/(sig-1);
gc2=ss_g/(sig-1);
gc3=ss_g/(sig-1);
gp1=-ss_g/(sig-1);
gp2=-ss_g/(sig-1);
gp3=-ss_g/(sig-1);
ga11=ss_g;
ga12=ss_g;
ga13=ss_g;
ga21=ss_g;
ga22=ss_g;
ga23=ss_g;
ga31=ss_g;
ga32=ss_g;
ga33=ss_g;
Ha11=log(ss_ha(1,1));
Ha12=log(ss_ha(1,2));
Ha13=log(ss_ha(1,3));
Ha21=log(ss_ha(2,1));
Ha22=log(ss_ha(2,2));
Ha23=log(ss_ha(2,3));
Ha31=log(ss_ha(3,1));
Ha32=log(ss_ha(3,2));
Ha33=log(ss_ha(3,3));
rp11=log(ss_rp(1,1));
rp12=log(ss_rp(1,2));
rp13=log(ss_rp(1,3));
rp21=log(ss_rp(2,1));
rp22=log(ss_rp(2,2));
rp23=log(ss_rp(2,3));
rp31=log(ss_rp(3,1));
rp32=log(ss_rp(3,2));
rp33=log(ss_rp(3,3));
eps11=log(ss_eps(1,1));
eps12=log(ss_eps(1,2));
eps13=log(ss_eps(1,3));
eps21=log(ss_eps(2,1));
eps22=log(ss_eps(2,2));
eps23=log(ss_eps(2,3));
eps31=log(ss_eps(3,1));
eps32=log(ss_eps(3,2));
eps33=log(ss_eps(3,3));
V11=log(ss_v(1,1));
V12=log(ss_v(1,2));
V13=log(ss_v(1,3));
V21=log(ss_v(2,1));
V22=log(ss_v(2,2));
V23=log(ss_v(2,3));
V31=log(ss_v(3,1));
V32=log(ss_v(3,2));
V33=log(ss_v(3,3));
J11=log(ss_j(1,1));
J12=log(ss_j(1,2));
J13=log(ss_j(1,3));
J21=log(ss_j(2,1));
J22=log(ss_j(2,2));
J23=log(ss_j(2,3));
J31=log(ss_j(3,1));
J32=log(ss_j(3,2));
J33=log(ss_j(3,3));
Vinnovator11=log(ss_vinnovator(1,1));
Vinnovator12=log(ss_vinnovator(1,2));
Vinnovator13=log(ss_vinnovator(1,3));
Vinnovator21=log(ss_vinnovator(2,1));
Vinnovator22=log(ss_vinnovator(2,2));
Vinnovator23=log(ss_vinnovator(2,3));
Vinnovator31=log(ss_vinnovator(3,1));
Vinnovator32=log(ss_vinnovator(3,2));
Vinnovator33=log(ss_vinnovator(3,3));
Jinnovator11=log(ss_jinnovator(1,1));
Jinnovator12=log(ss_jinnovator(1,2));
Jinnovator13=log(ss_jinnovator(1,3));
Jinnovator21=log(ss_jinnovator(2,1));
Jinnovator22=log(ss_jinnovator(2,2));
Jinnovator23=log(ss_jinnovator(2,3));
Jinnovator31=log(ss_jinnovator(3,1));
Jinnovator32=log(ss_jinnovator(3,2));
Jinnovator33=log(ss_jinnovator(3,3));
pi_share11=log(ss_PishareT(1,1));
pi_share12=log(ss_PishareT(1,2));
pi_share13=log(ss_PishareT(1,3));
pi_share21=log(ss_PishareT(2,1));
pi_share22=log(ss_PishareT(2,2));
pi_share23=log(ss_PishareT(2,3));
pi_share31=log(ss_PishareT(3,1));
pi_share32=log(ss_PishareT(3,2));
pi_share33=log(ss_PishareT(3,3));
pi_shareNT11=log(ss_PishareNT(1,1));
pi_shareNT12=log(ss_PishareNT(1,2));
pi_shareNT13=log(ss_PishareNT(1,3));
pi_shareNT21=log(ss_PishareNT(2,1));
pi_shareNT22=log(ss_PishareNT(2,2));
pi_shareNT23=log(ss_PishareNT(2,3));
pi_shareNT31=log(ss_PishareNT(3,1));
pi_shareNT32=log(ss_PishareNT(3,2));
pi_shareNT33=log(ss_PishareNT(3,3));
A11=log(ss_AT(1,1)*exp(T1));
A12=log(ss_AT(1,2)*exp(T1));
A13=log(ss_AT(1,3)*exp(T1));
A21=log(ss_AT(2,1)*exp(T2));
A22=log(ss_AT(2,2)*exp(T2));
A23=log(ss_AT(2,3)*exp(T2));
A31=log(ss_AT(3,1)*exp(T3));
A32=log(ss_AT(3,2)*exp(T3));
A33=log(ss_AT(3,3)*exp(T3));
Yw=ss_Yw;
end;
steady;
check(qz_zero_threshold=1e-20);
load Final_values.mat;
endval;
rho_chi1 = ss_rhochi(1);
rho_chi2 = ss_rhochi(2);
rho_chi3 = ss_rhochi(3);
B1=0;
B2=0;
B3=0;
distexo11=0;
distexo12=0;
distexo13=0;
distexo21=0;
distexo22=0;
distexo23=0;
distexo31=0;
distexo32=0;
distexo33=0;
distexoNT11=0;
distexoNT12=0;
distexoNT13=0;
distexoNT21=0;
distexoNT22=0;
distexoNT23=0;
distexoNT31=0;
distexoNT32=0;
distexoNT33=0;
chiexo11=0.;
chiexo12=0.0;
chiexo13=0.0;
chiexo21=0.;
chiexo22=0.0;
chiexo23=0.0;
chiexo31=0.;
chiexo32=0.0;
chiexo33=0.0;
r=log((1+ss_g/(sig-1))/(1+ss_g/(1-sig))/beta);
C1=log(ss_c(1));
C2=log(ss_c(2));
C3=log(ss_c(3));
Y1=log(ss_y(1));
Y2=log(ss_y(2));
Y3=log(ss_y(3));
YT1=log(ss_yT(1));
YT2=log(ss_yT(2));
YT3=log(ss_yT(3));
YNT1=log(ss_yNT(1));
YNT2=log(ss_yNT(2));
YNT3=log(ss_yNT(3));
T1=log(ss_TT(1)/ss_TT(3));
T2=log(ss_TT(2)/ss_TT(3));
T3=log(ss_TT(3)/ss_TT(3));
TNT1=log(ss_TNT(1)/ss_TNT(3));
TNT2=log(ss_TNT(2)/ss_TNT(3));
TNT3=log(ss_TNT(3)/ss_TNT(3));
P1=log(ss_p(1));
P2=log(ss_p(2));
P3=log(ss_p(3));
PT1=log(ss_pT(1));
PT2=log(ss_pT(2));
PT3=log(ss_pT(3));
PNT1=log(ss_pNT(1));
PNT2=log(ss_pNT(2));
PNT3=log(ss_pNT(3));
EX1=log(ss_EX(1));
EX2=log(ss_EX(2));
IM1=log(ss_IM(1));
IM2=log(ss_IM(2));
Pi1=log(ss_pi(1));
Pi2=log(ss_pi(2));
Pi3=log(ss_pi(3));
PiT1=log(ss_piT(1));
PiT2=log(ss_piT(2));
PiT3=log(ss_piT(3));
PiNT1=log(ss_piNT(1));
PiNT2=log(ss_piNT(2));
PiNT3=log(ss_piNT(3));
w1=log(ss_w(1));
w2=log(ss_w(2));
w3=log(ss_w(3));
Vinnov1=log(ss_vinnov(1));
Vinnov2=log(ss_vinnov(2));
Vinnov3=log(ss_vinnov(3));
chivar11 = log(ss_chi(1,1));
chivar12 = log(ss_chi(1,2));
chivar13 = log(ss_chi(1,3));
chivar21 = log(ss_chi(2,1));
chivar22 = log(ss_chi(2,2));
chivar23 = log(ss_chi(2,3));
chivar31 = log(ss_chi(3,1));
chivar32 = log(ss_chi(3,2));
chivar33 = log(ss_chi(3,3));
Z1=log(ss_Z(1));
Z2=log(ss_Z(2));
Z3=log(ss_Z(3));
Y_L1=log(ss_YL(1));
Y_L2=log(ss_YL(2));
Y_L3=log(ss_YL(3));
Hr1=log(ss_hr(1));
Hr2=log(ss_hr(2));
Hr3=log(ss_hr(3));
gz1=ss_g;
gz2=ss_g;
gz3=ss_g;
gt1=ss_g;
gt2=ss_g;
gt3=ss_g;
gtT1=ss_g;
gtT2=ss_g;
gtT3=ss_g;
gtNT1=ss_g;
gtNT2=ss_g;
gtNT3=ss_g;
gy1=ss_g/(sig-1);
gy2=ss_g/(sig-1);
gy3=ss_g/(sig-1);
gc1=ss_g/(sig-1);
gc2=ss_g/(sig-1);
gc3=ss_g/(sig-1);
gp1=-ss_g/(sig-1);
gp2=-ss_g/(sig-1);
gp3=-ss_g/(sig-1);
ga11=ss_g;
ga12=ss_g;
ga13=ss_g;
ga21=ss_g;
ga22=ss_g;
ga23=ss_g;
ga31=ss_g;
ga32=ss_g;
ga33=ss_g;
Ha11=log(ss_ha(1,1));
Ha12=log(ss_ha(1,2));
Ha13=log(ss_ha(1,3));
Ha21=log(ss_ha(2,1));
Ha22=log(ss_ha(2,2));
Ha23=log(ss_ha(2,3));
Ha31=log(ss_ha(3,1));
Ha32=log(ss_ha(3,2));
Ha33=log(ss_ha(3,3));
rp11=log(ss_rp(1,1));
rp12=log(ss_rp(1,2));
rp13=log(ss_rp(1,3));
rp21=log(ss_rp(2,1));
rp22=log(ss_rp(2,2));
rp23=log(ss_rp(2,3));
rp31=log(ss_rp(3,1));
rp32=log(ss_rp(3,2));
rp33=log(ss_rp(3,3));
eps11=log(ss_eps(1,1));
eps12=log(ss_eps(1,2));
eps13=log(ss_eps(1,3));
eps21=log(ss_eps(2,1));
eps22=log(ss_eps(2,2));
eps23=log(ss_eps(2,3));
eps31=log(ss_eps(3,1));
eps32=log(ss_eps(3,2));
eps33=log(ss_eps(3,3));
V11=log(ss_v(1,1));
V12=log(ss_v(1,2));
V13=log(ss_v(1,3));
V21=log(ss_v(2,1));
V22=log(ss_v(2,2));
V23=log(ss_v(2,3));
V31=log(ss_v(3,1));
V32=log(ss_v(3,2));
V33=log(ss_v(3,3));
J11=log(ss_j(1,1));
J12=log(ss_j(1,2));
J13=log(ss_j(1,3));
J21=log(ss_j(2,1));
J22=log(ss_j(2,2));
J23=log(ss_j(2,3));
J31=log(ss_j(3,1));
J32=log(ss_j(3,2));
J33=log(ss_j(3,3));
Vinnovator11=log(ss_vinnovator(1,1));
Vinnovator12=log(ss_vinnovator(1,2));
Vinnovator13=log(ss_vinnovator(1,3));
Vinnovator21=log(ss_vinnovator(2,1));
Vinnovator22=log(ss_vinnovator(2,2));
Vinnovator23=log(ss_vinnovator(2,3));
Vinnovator31=log(ss_vinnovator(3,1));
Vinnovator32=log(ss_vinnovator(3,2));
Vinnovator33=log(ss_vinnovator(3,3));
Jinnovator11=log(ss_jinnovator(1,1));
Jinnovator12=log(ss_jinnovator(1,2));
Jinnovator13=log(ss_jinnovator(1,3));
Jinnovator21=log(ss_jinnovator(2,1));
Jinnovator22=log(ss_jinnovator(2,2));
Jinnovator23=log(ss_jinnovator(2,3));
Jinnovator31=log(ss_jinnovator(3,1));
Jinnovator32=log(ss_jinnovator(3,2));
Jinnovator33=log(ss_jinnovator(3,3));
pi_share11=log(ss_PishareT(1,1));
pi_share12=log(ss_PishareT(1,2));
pi_share13=log(ss_PishareT(1,3));
pi_share21=log(ss_PishareT(2,1));
pi_share22=log(ss_PishareT(2,2));
pi_share23=log(ss_PishareT(2,3));
pi_share31=log(ss_PishareT(3,1));
pi_share32=log(ss_PishareT(3,2));
pi_share33=log(ss_PishareT(3,3));
pi_shareNT11=log(ss_PishareNT(1,1));
pi_shareNT12=log(ss_PishareNT(1,2));
pi_shareNT13=log(ss_PishareNT(1,3));
pi_shareNT21=log(ss_PishareNT(2,1));
pi_shareNT22=log(ss_PishareNT(2,2));
pi_shareNT23=log(ss_PishareNT(2,3));
pi_shareNT31=log(ss_PishareNT(3,1));
pi_shareNT32=log(ss_PishareNT(3,2));
pi_shareNT33=log(ss_PishareNT(3,3));
A11=log(ss_AT(1,1)*exp(T1));
A12=log(ss_AT(1,2)*exp(T1));
A13=log(ss_AT(1,3)*exp(T1));
A21=log(ss_AT(2,1)*exp(T2));
A22=log(ss_AT(2,2)*exp(T2));
A23=log(ss_AT(2,3)*exp(T2));
A31=log(ss_AT(3,1)*exp(T3));
A32=log(ss_AT(3,2)*exp(T3));
A33=log(ss_AT(3,3)*exp(T3));
Yw=ss_Yw;
end;
steady;
perfect_foresight_setup(periods=1000);
perfect_foresight_solver;
shock13 = [0:0.1:1];
shockchi31= log([1:0.1:1/0.4]);//log([1:0.1:1.6]);
shockchi32= log(1);//log([1:0.2:1/0.4]);//log([1:0.1:1.6]);
shockchi33= log([1:0.1:1/0.4]);//log([1:0.1:1.6]);
iterations = length(shock13)*length(shockchi31)*length(shockchi32)*length(shockchi33);
resultsC1 = zeros(iterations,1002);
resultsP1 = zeros(iterations,1002);
resultsgc1 = zeros(iterations,1002);
resultsC2 = zeros(iterations,1002);
resultsP2 = zeros(iterations,1002);
resultsgc2 = zeros(iterations,1002);
resultsC3 = zeros(iterations,1002);
resultsP3 = zeros(iterations,1002);
resultsgc3 = zeros(iterations,1002);
numit = 0;
for i=1:length(shock13)
    for j=1:length(shockchi31)
 for jj=1:length(shockchi32) 
for k=1:length(shockchi33)
                numit = numit + 1;                
                oo_.exo_simul(2:1002,7) = shock13(i);
                oo_.exo_simul(2:1002,21) = shockchi31(j);
                oo_.exo_simul(2:1002,24) = shockchi32(jj);
                oo_.exo_simul(2:1002,27) = shockchi33(k);
             
                oo_ = perfect_foresight_solver_core(M_,options_,oo_);
                if oo_.deterministic_simulation.status
                    resultsC1(numit,:) = oo_.endo_simul(strmatch('C1',M_.endo_names,'exact'),1:1002);//oo_.endo_simul(7,:);
                    resultsP1(numit,:) = oo_.endo_simul(strmatch('P1',M_.endo_names,'exact'),1:1002);//oo_.endo_simul(14,:);
                    resultsgc1(numit,:) = oo_.endo_simul(strmatch('gc1',M_.endo_names,'exact'),1:1002);//oo_.endo_simul(27,:);
                    resultsC2(numit,:) = oo_.endo_simul(strmatch('C2',M_.endo_names,'exact'),1:1002);//oo_.endo_simul(34,:);
                    resultsP2(numit,:) = oo_.endo_simul(strmatch('P2',M_.endo_names,'exact'),1:1002);//oo_.endo_simul(41,:);
                    resultsgc2(numit,:) = oo_.endo_simul(strmatch('gc2',M_.endo_names,'exact'),1:1002);//oo_.endo_simul(54,:);
                    resultsC3(numit,:) = oo_.endo_simul(strmatch('C3',M_.endo_names,'exact'),1:1002);//oo_.endo_simul(61,:);
                    resultsP3(numit,:) = oo_.endo_simul(strmatch('P3',M_.endo_names,'exact'),1:1002);//oo_.endo_simul(68,:);
                    resultsgc3(numit,:) = oo_.endo_simul(strmatch('gc3',M_.endo_names,'exact'),1:1002);//oo_.endo_simul(81,:);
                    position(numit,:)=[numit i j jj k];
                end
end
end
            end
        end
